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The aim of this lecture is to study numerical algorithms for the integration of 
stochastic differential equations. I will derive an algorithm which exactly integrates 
the SDE, using a generalization of a Taylor series in the presence of stochastic 
forces. Given the complexity, we will find that this algorithm, although "exact", it 
is not particularly fast or, in general, the most convenient in a digital simulation, 
although it is a useful benchmark to test other algorithms, and I will try to improve 
it. I will discuss the features of different algorithms, both in terms of accuracy in a 
deterministic sense, and also in statistical terms, i.e. how well the algorithm is able 
to reproduce, for instance, the correct equilibrium distribution. I will then briefly 
introduce algorithms to integrate stochastic differential equations which are driven 
by correlated noise. The lecture will close with the discussion of a few algorithms 
for the special case of a two dimensional system in a potential, subject to damping 
and noise. The interest for this particular case is due to the importance and the 
interest in algorithms which can integrate, for instance, particles in the liquid state. 

1 Basic algorithm and relation to Fokker-Planck equations 

1.1 Introduction 

A very simple and straightforward algorithm to numerically integrate one di- 
mensional differential equations driven by a single stochastic force (external 
to the system) was introduced by Rao in [1]. Note that some of the technical 
terms used in this introduction will become clear in the rest of the lecture. 
The algorithm is really a Taylor expansion, or, more precisely, a one step col- 
location scheme. It is possible to use other integration schemes (for a review 
of some widely used algorithms see Mannella in [2]; more recent papers are 
references [3-11]). To integrate stochastic differential equations (SDE), we 
can basically have integration schemes based on predictor correctors [12-19] 
and schemes based on Runge-Kutta [20-25,4,5]. Interesting material can also 
be found in the papers by Riimelin [26] and Riggs Jr. [27,28] where, beside 
Runge-Kutta based approaches, other integration methods are presented and 
some some attention is devoted to the accuracy of the different schemes. Fi- 
nally, other particular schemes, not immediately connected to the two classes 
just introduced are the ones proposed in [10,11]. Note that I will only deal with 
temporal noise: for SDE driven by temporal and spatial noise, the interested 
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reader could refer to, for instance, [29] and references therein. 

In the literature it is possible to find algorithms which share the same 
basic idea with the one I will discuss here: see for instance [30-36]. Yet a 
different approach is the one followed by Fox [37], where the deterministic force 
appearing in the sde is integrated via a Runge-Kutta whereas the stochastic 
force is integrated using a Taylor expansion The algorithms proposed in 
references [4,5] extend the Fox approach to an implicit Runge-Kutta treatment 
of also the stochastic force. 

That are some reasons which make one step collocation schemes slightly 
preferable to predictor correctors and Runge-Kutta approaches. Predictor cor- 
rectors of high orders, which by definition imply continuous and differentiable 
functions, seem to be of difficult or risky application to SDE which, by def- 
inition, have continuous but non-differentiable solutions. In practice, the al- 
gorithm one derives yields a trajectory which is indeed non-differentiable, but 
the derivation of the schemes seem to be rather empirical. We should add here 
that at low order, when it is possible to make a comparison, typical predic- 
tor correctors coincide with schemes based on Taylor expansions. As far as 
Runge-Kutta are concerned, we have similar problems to the ones mentioned 
for predictor correctors, with the added complication that schemes are known 
in the literature which in principle should coincide but in practice differ, the 
differences being due to different ways of expanding the equilibrium distribu- 
tion of the stochastic force. There is however a more fundamental reason to 
prefer one-step collocation schemes. In this approach, for relatively small or- 
ders of the schemes, it is possible to simply evaluate the different stochastic 
integrals which appear in the SDE, with the bonus that the algorithm can be 
straightforwardly generalized to cover cases for which the stochastic force is 
non-white. 

I will always use the Stratonovich calculus [40] to integrate the stochastic 
integrals which will be needed at the various stages. After the derivation 
of the algorithm, a section will be devoted to explain the details of using a 
given calculus and the relation between a stochastic differential equation and 
a Fokker-Planck equation. 

1.2 Derivation of the algorithm 

The most general SDE in which the stochastic force is linear can be written in 
the form 

ii{t) = fi{xi{t),t)+gi{xi{t),t)C{t). (1) 

^According to private communications from R.F. Fox, the Runge-Kutta is normally a 
fourth order RK [38,39]. 
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I will always assume that ^(t) is a Gaussian variable. To describe it, then, 
I will need only the first two moments. For simplicity, at least in this section, 
these moments will always be given by (the symbol (...) will denote averages 
over the noise realizations) 

m) = (2a) 

mm) = Ht-t'), (2b) 

where S{t) is the usual Dirac delta. Note that Eq. 2a also defines the spectral 
distribution of the stochastic force: the spectral distribution of the stochastic 
force is given by the Fourier trasform of the two-time correlation function. For 
this case, the two-time correlation function is a Dirac delta, hence the spectral 
distribution will be trivially fiat (the so called white noise). The derivation of 
the basic algorithm is much more simple in this case. For the moment, I will 
assume that and gi are autonomous, i.e. they are not explicit functions of 
time. 

Expand via the Taylor formula the functions in Eq. 1. It is possible to 
write (I will assume a summation over repeated indices) 

Mxj(t)) = Mx,(0)) + djMxm)(xj(t) - xM) + ■■■ (3a) 
giixjit)) = giixjiO)) + djgiixkiO)){xj{t) -Xj{0)) + ..., (3b) 

plus derivatives of higher orders. I can also rewrite the above equation as 

f' = f° + f°J^Ai) + lf°,,J^Ai)^Mi) + ■■■ (4a) 
al = gi+gijSxj{t) + ]^g°j,,Sxj{t)Sxk{t) 

+ ^^9ijki^Xj{t)Sxk{t)Sxi{t) + ..., (4b) 

where, for instance, g^ jf; = q^.q^^^ 9i(xi (i = s)) and so on: the symbol 5xi(i) 
means Xi[i) — Xi[Q) . 

Eq. 1 can be integrated formally writing 

Xi{t) - Xi{Q) = i Xi{s)ds = i fi{xi{s))ds + / gi{xi{s))^{s)ds. (5) 
Jo Jo Jo 

The integration is formal because the rhs of Eq. 5 still depends explicitly on 

Xi{t). 

The idea is now to integrate Eq. 5 for times between zero and h, where h 
is a (small) integration time step, then to substitute the Taylor expansions of 
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/* and gj in the r.h.s. of equation 5. Treating h as an expansion parameter, 
it will be possibile to derive higher orders approximations for Xi{h) — Xi{0), 
approximations which will be substituted back in the Taylor expansions of 
/* and (jf*, generating terms of higher orders for Xi{h) — Xi{0) and so on. A 
symbol widely used in the following will be Sxf{h) indicating the contribution 
to Xi{h) — Xi{0) coming from the perturbation term of order h" . 
It is now possible to write 

Sx,{t) = ^ {f°+ f°Jx,{s) + lf°,Jx,{s)Sxk{s) + . . .)ds 

+ £,{s){9° + glj5xj{s) + -g°jk5xj{s)5xk{s) 

+ Y\^ljkiSxj{s)5xk{s)5xi{s) + . . .)ds. (6) 

At the lowest possible order (i.e. h^^^), obtained keeping only the term g'^ 
in the second integral appearing in equation 6, one has straightforwardly 

Sxy\h)=g°Z,ih)+oih^) (7) 
Zi{h) = / C{s)ds = VhYi (8) 



where Yi is a stochastic Gaussian variable with average zero and standard 
deviation one; Eq. 8 follows considering that Zi{h) = J £,{t)ds is a linear 
combination of gaussian random variables, hence a gaussian variable. The 
moments of Zi{h) follow, having used Eq. 2a to compute the averages. It is 
very simple to show that x^^^h) in Eq. 7 is of order h^^^: it follows from the 
definition of Zi{h). 

Now, let me substitute equation 7 into equation 6, keeping only the lowest 
order terms (i.e. 0{h^)): 



Sxj(h)= l' f°ds+ f\(s)gl,Sxl^\s)ds = f°h + gl,g° f Z,{s)i{s)ds 
Jo Jo Jo 

f°h + lglkgl[Zi(h)? + o(h'^'), (9) 



where the equality 



Z,is)as)ds=lz,ihf 

^ 
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follows from the definition of Z\ (h) (Eq. 8) and from the Stratonovich calculus 
for the stochastic integrals (see below). 

The details of how the calculations follows can be found in Mannella [2], 
including a derivation of the various stochastic integrals which are needed, i.e. 

Z,{h) ^ £ Z,{s)ds = /^t 1^ + ^1 (10) 

and 




where Y2 and Y3 are two uncorrelated gaussian deviates with average zero and 
standard deviation one. 

The Taylor expansion around t = Q reads (see Mannella in [2] for the 
intermediate algebra) 

x,(h) = x,(0)+g°Z^(h) +f°h+ lglkgk[Zi(h)f + Z2(h) {fl,g° - glJ°} 
+ hZ,{h)glJ° + l^[Z,{h)f {gl,g°g°, + gl,gl,g°} 
+ \kjkg°9lZ:,{h) + gl, {fl^g° - gljf} {Z,{h)Z,{h) - Zs{h)} 

+ Iflk {fkh'+gl,gp3(h)} + l^gl, {giyr.,9° + 9l,n9°9n} [Zi(h)T 

+ 9lAj: [\[Zi{h)f - 2^3(/»)} + \9l,,9',f, {h[Zi{h)f - Z:,{h)] 
+ Y^9l,,9ln9'n9'AZi{^)f + l9l,,n9',9Wn[Zi{h)]'' + o{h^l^) (12) 

which is the final expression. It will be will shown further down in some 
selected examples that the algorithm of Eq. 12 is good enough for a general 
purpose integration, and that higher orders in h are not really needed. On 
the other hand, to compute terms at (even) next perturbation order would be 
a very difficult task, due to the appearance of non-Gaussian stochastic terms 
(for instance, like Z3 above) for which the appropriate statistics can only be 
approximated. 

It is very simple to generalize the basic algorithm to take into account 
non-autonomous functions. It is enough to add to Eqs. 4 the terms due to 
the partial derivatives with respect to time of and of gi, i.e. the terms 
fi^) fli^ and fli j{xj{t) — Xj{0))t (higher order terms will not contribute at the 
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perturbation order considered in Eq. 12). In practice 



Sx^^\hy =g°{hZ,{h)-Z2{h)) (13a) 

S^lihY = y/° + Ig^l {h[Z,{h)f - Zsih)} (13b) 

+ l9lk9k {h[Zi{h)f + Zs{h) - 2Z,{h)Z,{h)} (13c) 

are all the necessary corrections to Eq. 12 of order 0{h^). 

Note that the algorithm was derived under the assumption of only one 
stochastic forcing: the case of a multidimensional stochastic forcing is much 
more complex and we refer the reader to Mannella in [2] for further details. 

1.3 Fokker-Planck equations, the Ito Stratonovtch controversy and numerical 
simulations 

As strange as it may sound, the Ito Stratonovich controversy is anything but 
a controversy [41]. For some time, within the stochastic physics community, 
there has been a debate on which should be the correct prescription to obtain 
the Fokker-Planck operator corresponding to a given Langevin equation. Let 
me clarify the problem with an example. 
Suppose I have the Langevin equation 

x = f{x)+g{x)^{t) (14) 

where ^(t) is as usual a Gaussian variable with moments 

m) = (15) 

mm) = ^DS{t - s) (16) 

it is in general possible to write an infinite number of different Fokker-Planck 
equations (the differential operator driving the probability distribution of the 
variable x). The most common prescriptions found in the literature are that 
due to Ito [42], and that due to Stratonovich [40]. In particular, if P{x,t) is 
the probability distribution of x as function oft under the fiux given by Eq. 14, 
we obtain (I (S) will refer to Ito (Stratonovich)) 

l^Pi^,i) = ^ {-/W + l9ix)g'{x) + Pi^,i) (S). (18) 
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How can it be possible that given the same differential equation (Eq. 14) we 
should have two (better, infinite, as we will see) different differential operators 
for the evolution of the probability distribution, if (note the structure of Eq. 14) 
g{x) depends explicitly on x7 The solution is that Eq. 14 does not make sense 
per se. It is a fact that the stochastic process ^(t) in Eq. 14 is such that the 
quantity x is "kicked" by the noise, with the result that its integral {x{t)) is not 
continuous: in practice, it is impossible to evaluate the term g{x)£,{t) because 
it is unknown for which x g{x) should be computed. 

The problem has to do with the intrinsic limits of the Riemann-Stieljes 
integral, in the presence of stochastic forces [43,44]. Let me suppose that I 
have two functions f{t) and g{t), defined in the interval a < t < h. For each 
partition P : a < to < ti < . . . < t„ = h I can build 

n 

Sp = J2fi^')i9{ti)-g{U-i)], (19) 

i = l 

where < ^j- < ti. Define |-P| = maxi<j<„(tj' — tj_i). If it exists, and it is 
finite, the limit 

lim Sp = S, (20) 
|PKo 

I will term S as the Stteltjes integral of f[t) with respect to g{t). I will sym- 
bolically write it as 

S = I fit) dg{t). (21) 

Now, if I have [44] that g{t) is the difference between two functions which are 
finite, monotonic and non-decreasing, and also that f{t) is continuous, then it 
follows that the integral S exists. Also, if J^* f{t) dg{t) exists, then the integral 
fa df{t) exists too, and in particular I have 

fit) dg{t) = f{b)g{b) - f{a)g{a) - f g{t) df{t). (22) 

Finally, ii g[t) is differentiable and g' {t) and f[t) are integrable, I have also 
that ^ ^ 

' f{t) dg{t) = [ fit) g'{t) dt. (23) 

Let me now introduce a new quantity, Wt, known as Wiener process. Wt 
is defined (see Eq. 14) as 

Wt= [ C{s) ds, (24) 
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or, symbolically but formally less accurately, 



dWt = £,{t) dt. (25) 

Eq. 23 implies that if I had to evaluate the integral g[s) ^(s) ds, I could 
instead evaluate the integral 

g{s)i{s)ds= / g{s)dW, = g(t)Wt - / W,g'{s)ds (26) 

JO Jo 

which obviously is defined only if g{s) does not depend on Wt, otherwise the 
function g' {t) is not defined. 

The solution of this problem is to extend the original definition of Riemann- 
Stieltjes integral to include the case of stochastic functions. Let me first intro- 
duce the characteristic function of the interval [a, &], which I will call X[afi]{'t), 
defined as 

, , J 1 if a < t < & , , 

X[aMW = [ otherwise (^7) 

Let f[i) be a step function in the interval [a,&], i.e., if I have the partition 
a <tQ <ti < . . . <tn = i, 

n-1 

/W = E/(^Ox[t.,t.+.]W- (28) 

i = 

I will define (Ito) 

»6 n-1 

(I) / f{t)dWt = Y.f{ti)[Wt,^,-Wt,]. (29) 

i=o 

Eq. 29 is the definition of integration according to the Ito prescription. We 
note that the function f[t) could in principle be a stochastic function, given 
that in the Ito prescription only the values it assumes at the times ti matter. 
Physically, we could interpret Eq. 29 saying that to evaluate the stochastic 
integral I must take the value the function f[t) assumes mimedmtely before the 
stochastic term is "applied" , to overcome the problems with the discontinuity 
of the variable t. From Eq. 29 I can obtain the following identity [43,44] 
(integration by parts) 

/6 1 ^ oo ^ ^ 

WsdW. = - {Wi - Wl) -TY.^dWnf = j iyVl - Wl) -j[h-a). 
n = 

(30) 
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It should be clear that the definition introduced in Eq. 29 is somehow 
arbitrary. Strictly from a mathematical point of view, the Ito integral has 
interesting properties (see [43]) which makes it formally most attractive. From 
a physical point of view, on the other hand it is clear that Eq. 29 is not 
symmetric with respect to the variable t (somehow it "points towards the 
future"). Also, Eq. 30 show that the usual rules for parts integration do not 
apply to the Ito definition. A different prescription, which is symmetric with 
respect to the variable T and for which the usual part integration rules apply 
is the one introduced by Stratonovich [40]. Using the same hypothesis used for 
Eq. 29 we can define 



»6 n-l 

(S) / m dW,^J2 ''^"\''^'^" (31) 
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Clearly, Eq. 31 is symmetric with respect to the variable t; also, we have that 



(S) 



(32) 



In general, it is even possible to think of some "mixed" prescription, which 
could be defined as 



»6 

(M) / f{t)dWt = Y, 

•'a -n 



M 



+1J 



- Wt,] , (33) 



for which we have 



(M) 



dW, 



2-^ 



(&■ 



(34) 



which yields the Ito (Stratonovich) prescription taking M = {\) . 

Typically a numerical algorithm (and in particular the one step collocation 
derived here, as mentioned before) uses heavily a standard rule of integration by 
parts for the evaluation of the stochastic integrals at different orders. This im- 
plies that the algorithm integrates according to the Stratonovich prescription. 
If it were necessary to integrate according to the M prescription (remembering 
that for Ito we have M = 0), it would be enough to replace in Eq. 1 with 

ji + (M — \) (^gf^fi'i^ fi'j) assuming a summation of repeated indices. Let me 
stress that a numerical simulation cannot, under any circumstance, resolve the 
Ito Stratonovich controversy, because the actual form given to the integrating 
algorithm must be decided beforehand: and if a prescription M is picked to 
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write the algorithm, obviously the integration scheme will integrate according 



It is possible to say that the Stratonovich prescription corresponds to con- 
sidering a stochastic process with a finite correlation time, where the limit of 
zero correlation time is taken once all the different integrals have been com- 
puted. This poses the following problem: in practice, what is the "correct" 
prescription to integrate Eq. 14? A possible answer is the following: if the 
"internal" fiuctuations of the variable x are faster than the fiuctuations of the 
variable ■J(t), then it is reasonable that the correct prescription will be more Ito 
like. In the opposite limit, the Stratonovich prescription will be the norm. This 
has been shown to be the correct picture in some analogue simulations [45], 
where, electronically simulating the Langevin equation 



(where ^(t) and are two uncorrelated stochastic processes) it is possible 
to show that the equilibrium distribution of the variable x will move from the 
one characteristic of Ito prescription to the one typical of the Stratonovich 
prescription simply changing the relative (short) correlation times of the two 
stochastic processes. 

1.4 Improving the basic algorithm 

In principle it is possible to improve the basic algorithm combining it with a 
suitable number of predictor-correctors [38]. The first attempt in this direction, 
following the approach we have shown for the derivation of the stochastic 
variables in the basic algorithm, it is probably the one by Blum [46]. He 
introduced the so called Heun scheme, which has been subsequently used also 
by Riimelin [26]. In practice the method consists of an Adams Bashforth (AB) 
predictor of order zero corrected via an Adams Moulton (AM) corrector of 
order one, where the basic algorithm at first order is used to evaluate the 
forces in both stages (see Eqs. 7 and 9). We would like to point out that the 
introduction of an AM corrector, which weights equally the stochastic forcing 
at the beginning and at the end of the integration time step, will automatically 
implement the Stratonovich calculus for the SDE (see previous chapter for more 
details). It would not be necessary, then, to introduce the term jgj[Zi{h)]''^ 
(see Eq. 9), which in the basic algorithm made sure that the calculus employed 
was really the Stratonovich one. 

Let me define AB(n) (AM(n)) as an Adams Bashforth predictor (Adams 
Moulton corrector) of order n. More in details, having to integrate the SDE 



to M. 



X = 



f{x)+g{x)^{s) + 7j{t) 



(35) 



X = 



f{x)+g{x)^ 



(36) 
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m) = (37) 

mm = Ht-s), (38) 

with ^ Gaussian, according to the Heun algorithm at each integration step one 
will have to compute 

1. £{h) = x{0) + f{x{0)) + g{x{0))Zi{h) 

2. x{h) = x{0) + ^ [f{x{0)) + f{£{h))] + ^ [g{x{0)) + g{£{h))] Z,{h).{39) 

and go through the loop again for the following step of integration. 

Independently, Mannella and Palleschi [47] used an algorithm which con- 
sists of an AB predictor of zero-th order, corrected via an AM corrector still 
of zero-th order, but using, as the building block for the forces, the basic al- 
gorithm at order h^^^. In this case, however, it is important to note that the 
corrector requires the evaluation of the stochastic integral corresponding to 
Z2{h) around h. The corresponding quantity is evaluated in [48]. However, 
this algorithm is here quoted only for the sake of completeness, because from 
the point of view of reconstructing the correct equilibrium distribution it does 
not do better than simpler and faster algorithms. 

2 Behavior of the different integration schemes 

As mentioned in the previous pages, when considering the integration of SDE, 
one has two different aspects: the accuracy of the integration of the "determin- 
istic" part of the SDE and how well represented are the statistical quantities 
associated with the given SDE. 

In the following I will use two tools to study the different schemes. I will 
use, for the deterministic drift, the standard tool normally employed to work 
out the coefficients in a predictor corrector scheme (it will be clear in the fol- 
lowing how the technique works). For the statistical properties, I will work 
out the actual equilibrium distribution obtained in the numerical integration: 
and from a comparison with the exact quantities we will easily judge pros and 
cons of the various schemes. This is achieved deriving the Fokker-Planck equa- 
tion obtained in the numerical scheme. Suppose that the discrete integration 
scheme reads 

Xi{t + h) = Xi{t)+Fi{x,t) (40) 

(see the structure of Eq. 12), then the following partial differential equation 
follows (P is the probability distribution) 

Pii+k)-Pii)=j: j: ^...^AY^.i^w (41) 

n = l Xi...Xn 
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where 



nl 



(42) 



This equation is obtained [49] from the evolution equation 



P{xi,t + h) = { / dyi . . .dyr,P{yi,t)Uy^S{xi{t + h)-yi{t)-Fi{y,t)))r,oise (43) 



which means that we can divide both members of Eq. 41 by h. In the limit 
/i — 7- the l.h.s. of Eq. 41 will give the derivative of P{t) with respect to 
time. However, depending on how the various powers of h will vanish in the 
r.h.s. of Eq. 41 , we will infer which integration scheme is preferable at finite 
integration time steps. 

I will look at the following cases: 

Heun scheme: The Heun scheme was described in Eq. 39 

Simple Euler: The simple Euler corresponds to the algorithm given by 
Eqs. 7 and 9. 

Exact propagator: I will define as the exact propagator algorithm an algo- 
rithm where I integrate exactly (i.e., with a very high order integration 
time step) the determimsUc part of the equations of motion, and then 
the contribution given by Eq. 7 is added to the deterministic fiow. 

Algorithm o[h'^): It will be the integration scheme given by Eq. 12, dropping 
the 0{h^) terms. 

Zs-less: It will be the algorithm given by Eq. 12, but without the terms 
containing Z^. 

Full algorithm: The algorithm given by Eq. 12. I will also refer to this case 
as to the algorithm 0{h^). 

We are now ready to compare the different algorithms. 



by Taylor expansion and average over the noise realization. 
The r.h.s. in Eq. 41 will turn out to be always of the form 



oo 




(44) 



n = l 
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2.1 DetermimsUc behavior of the different integration schemes 



The evaluation of the error associated with a given integration scheme parallels 
the derivation of a predictor corrector scheme. To understand how the thing is 
done, I will derive the AM corrector to first order. Starting from the equation 

x = f{x), (45) 

the first order AM reads 

x{t + h) = x{t) + h{aif{x{t + h)) + aof{x{t))). (46) 

We need to evaluate (to find the optimal values) for the parameters ao and ai. 
The idea is to write 

x(t + nh) = (t + nh)"" (47) 

where m is by one unity larger than the unknowns one has to optimize: in this 
case, m = 3. 

Inserting Eq. 47 in Eq. 46, noting that f{x{t + nh)) = x{t + nh), I have 

{t + nh)^ =t^ + h{3ai{t + h)^ + 3aot^) (48) 

from which it follows 

t^ + it^h + ith'^ + h^ = t^ + t^h{iai + 3ao) + Uh'^ai + ih^ai. (49) 

Equating now the terms which have the same power oft, we immediately find 
the standard result, ai = a^ = 1/2. However, it is evident that after all 
algebra, we are left not with zero, but with the quantity 

/i3 = 3/2/1^ (50) 

which implies that the intrinsic error associated with this method will be 
0(/iV2). 

For a simple minded (Euler) integration one has 

(t + /i)2 = ^2 _^ 2M (51) 

which leads to an intrinsic error which is 0{h^). 

The algorithm of Eq. 12, if restricted to the deterministic part, implies 

that 

xit + h) = xit) + hfixit)) + y/(*W)^^^^- (52) 
13 



Noting that x = f{x{t)) and that x = xf'{x{t)) = f{x{t))f'{x{t)), it follows 
that in this case 

x{t + h) = x{t) + hx{t) + h'^x/2 (53) 
and, using x(t) = , I have 

(t + hf = + m'-' + 3hH (54) 

which leads to an intrinsic error which is 0{h^). This error is of the same order 
of the AM corrector I examined before, although it is somehow larger (note 
that in the AM corrector example we had a factor 1/2 in front of h^). We 
should expect, hence, that the AM corrector is slightly more accurate than the 
algorithm of Eq. 12 in the limit of very small noise intensities. 

I will now verify these results in a simple stochastic system. The system 
considered in the numerical simulations is a model of dymer in the presence of 
thermal fluctuations [50-52]: it is the semiclassical approximation of a nonlin- 
ear Schrodinger equation which can be cast in the form of a classical Langevin 
equation. This dymer, incidentally, is basically equivalent to a spin in the pres- 
ence of a magnetic held along the z axis and of a fluctuating magnetic fleld on 
the X axis. The SDE one writes is in the form [51] 

i> = 2Vp (55a) 

q = —2Vp — x^P ~ X^^ (55b) 

r = xqp + xqz (55c) 

i = -Tz +C{t) -2Vq (55d) 

m) = 

where the different p, q are r are suitable linear combinations of the elements 
of the dymer density matrix, and the stochastic force ■J(t), which is assumed 
to have a Gaussian statistics, represents the interaction between the dymer 
and a thermal bath with flnite temperature. Although the last equation in 55 
is reminiscent of a colored noise evolution, the presence of the reaction term 
—2Vq in the equation for i implies that colored noise algorithms cannot be 
used. In practice it would be straightforward to modify these algorithms to 
cover this case, given the linearity of the equation for i, but here I really want 
to learn something about the different "deterministic" schemes, and I will then 
simply use standard white noise algorithms. 

Apart from the wide relevance of Eq. 55 in the physical sciences (the 
onset of possible localized states, for instance, could have consequences on the 
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onset of Davydov solitons in molecular chains), their importance here is that 
they admit an exact integral of motion: in fact, despite the nonlinearity, by 
inspection it is possible to note that the quantity p^+g^+r^ is an exact constant 
of motion (remembering that the equations are some kind of semiclassical 
approximation, it is possible to show that the constant equals one). For the 
numerical experiments presented further down (figure 1) the parameters chosen 
were x = 1-25, F = 3, 2V = 1, starting from the initial condition p = 1, q = 
r = 0. This implies that for the subsequent evolution + + = 1. A test 
particle was then followed for times equal to 200, and the final value of the 
quantity + + — 1 was plotted as function of the integration time step, 
for different integration algorithms and for two different noise intensities D. 




Figure 1: Algorithms for the dymer model, with D = 1.0 (left) and D = 0.01 (right). O 
exact propagator, + Simple Euler, □ Heun scheme, X Algorithm o{h?)^ A Algorithm 0{h?). 

For both intensities of the noise considered, the algorithms which are first 
order in the deterministic force are considerably less accurate in preserving the 
norm. Of the higher order algorithms, the Heun scheme does better for small 
noise intensities, whereas the full algorithm is better for large noise intensities. 
The algorithm which integrates exactly the deterministic force is better than 
the first order algorithms, but it is inferior to both the Heun scheme and the 
full algorithm. 

2.2 Statistical behavior of the different integration schemes 

The tool used to understand the statistical error associated with each integra- 
tion scheme is Eq. 41, where I will insert the numerical scheme in the expression 
for A'l „ and average over the noise variables. 

I will also limit my investigation, without loss of generality, to one dimen- 
sional system driven by one additive noise. The general system I will look at 
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is in the form 

x = -V'{x) + f{t) 
where f{t) is a gaussian random process, with moments 

{fit}) = {f{t)f{s)) = 2DS{t - s) 



(56) 



(57) 



Associated with the stochastic differential Eq. 56 we also have a Fokker- 
Planck equation, which reads 



dP{x,t) 
dt 



d_ 
dx 



V'{x) + D 



d_ 
dx 



PixA) 



Note that this equation is in the form of a conservation equation 

dP{x,t) dJ{x,t) _ ^ 
dt dx 

with the "probability current" J[x,i) given by the equation 

d 



J(x, t) 



V'{x) + D 



da 



P xj) 



(58) 



(59) 



(60) 



The equilibrium solution of Eq. 58 is readily obtained. First, given that 
we are talking of an equilibrium solution (I will also use mild hypotheses, like 
that the function V{x) grows quickly enough for large x's, so that a ground 
state for P[x,i) is defined), I can say that 



with the equation 



dP(x,t) 
dt 



must vanish. This leaves me 



d_ 
dx 



V'{x) + D 



PixA) 



0. 



(61) 



This equation implies that the divergence of the current must vanish: this 
implies that the current itself must be a constant. However, under the mild 
hypotheses mentioned above, the current must also vanish for \x\ oo, which 
implies that the current must vanish identically for any x once the equilibrium 
solution has been reached. 

This leads to the equation 

V'{x) + D^ P{x,t) = (62) 

which has the solution 



P{x,t) = Nexp{-V{x)/D} 



(63) 
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where N is a normalization constant. 

In general, integrating Eq. 56 with a discrete integration routine we will 
have an equilibrium solution which is not exactly given by Eq. 63, but rather 

by 

P{x,t) = Nexpi^(^-V{x) + J2h"S„{x)^ (64) 

In the figures which follow I have simulated, using the difli"erent schemes, 
the stochastic dynamics of the systems 

i = -x-x^ + f(t) (65) 

and 

i = x-x^ + f{t), (66) 

and I will then compare the result of the simulations with the "wrong" equi- 
librium distribution which is the equilibrium distribution obtained using each 
integration scheme with a finite time step. In the simulations, f{t) is a gaussian 
random variable as in Eq. 56. The comparison between theory and numerical 
simulations is shown in figure 2. 




0.0 -0.2 0.0 0.2 

X X 



Figure 2: Equilibrium distributions obtained using different integration scliemes, D = 
0.1, h = O.f. X and dotted line: exact propagator; A and dashed line, simple Euler; +, 
algorithm 0[h^); □, tleun algorithm; O? .^^3-less; Solid line, exact equilibrium distribution. 
Left, Eq. 66; right, distribution near the maximum for Eq. 65. 



2.3 Simple Euler 

In this case, the integration scheme reads 

x{t + h) = x{t) + hf{x{t)) + V2DhT] (67) 
17 



and we need only to keep the lowest four terms in Eq. 42. Carrying out the 
necessary algebra, the equilibrium distribution reads 




Note that the difference with respect to the true equilibrium distribution is 
order of h in the exponent. It is clear from the figures that indeed the simula- 
tions carried out using a simple Euler scheme follow very closely the expected 
theoretical distribution. 

2.4 Exact propagator 

In this case, I need to first solve the equation 



over the time interval [t,t + h], starting from x = x{t); call x{t + h)det the 
solution at t = t + h. Then, the stochastic x{t + h) is given by 



As above, it is possible to insert this prescription in the equation for the 
moments 42. After carrying out the necessary algebra, the equilibrium distri- 
bution reads 



which again differs from the true equilibrium distribution by a term order of h 
in the exponent. Incidentally, it is clear that at very small D the equilibrium 
distribution obtained using an exact propagator is worse than the equilibrium 
distribution obtained using a simple Euler scheme, as it is clear both from the 
figure and by inspecting the corresponding equilibrium distribution. 

This should already be a very important warning: when dealing with SDE, 
it is extremely delicate how higher orders are introduced; and if the orders in 
the approximation in the deterministic and in the stochastic part of the fiow are 
not "balanced" , when we increase the order of the algorithm we may actually 
make it worse! 



X = 



(69) 



x{t + h) = x{t h)det + V^Dhrj. 



(70) 
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2.5 Algonthm o(h'^) 



In this case we find a result similar to the one found for the simple Euler case: 
the equilibrium distribution generated by this scheme is given by 



P(xJ) 



N exp 



-V{x) + -{V'{xr/2-DV"{x)) 



IK 



(72) 



2.6 Zs-less 



In this case we find already an improvement with respect to the previous cases: 
the equilibrium distribution reads 



Pix,t) 



N exp 



-V{x) - 2 iDV"{x)) 



/D 



(73) 



In other words, the term proportional to V'{x)^ disappears. This implies 
that, comparing the results for the Euler scheme, the exact propagator, the 
algorithm o[h'^), and Eq. 12, the term V'{xY is canceled both by the joint 
"action" of the term Z2 and by the term proportional to fi^kfk in Eq. 12. But 
it also means that dropping one of these two ingredients yield a much less 
accurate integration scheme. Although I will show below that both the Heun 
scheme and the full algorithm are actually better integration scheme, for very 
small D the Z^-less algorithm is a significant improvement with respect to the 
other algorithms seen so far. 



2. 7 Heun scheme 

The algebra corresponding to this case is somehow more cumbersome, because 
now the correction to the true equilibrium distribution is order of h'^ in the 
exponent. Writing the equilibrium distribution in the form 

P{x, t) = N exp { [-V{x) + h'^F{x)] / D} , (74) 

the function F{x) satisfies the equation {V = V{x), F = F{x)) 

= UD'-'F' + 2V'^ - ADV'^V" + ?,?,D'^V'V"'^ + 21D'^V''^V"' 

1 - 24^3 y'V" - 2^3^ + 3D4y(^) . (75) 

The importance of the Heun scheme, however, is that the correction is only 
order of h'^ in the exponent, which means that the equilibrium distribution it 
generates will in general be fairly close to the true one, as it is confirmed by 
the simulations. 
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dh 



Figure 3: Left: MFPT with different algorithms, for D = 0.1, vs dh. O? full algorithm; 
□ , exact propagator; +, algorithm o[h^); X, Heun algorithm; A, simple Euler. Right: 
Comparison between the MFPT obtained using the white noise Heun algorithm (A), and 
the MFPT obtained using the coloured noise Heun algorithm, with r = 10~^ (O)- 

2.8 Full algorithm 

As in the case of the Heun scheme, the algebra in this case is very cumbersome. 
The correction to the true equilibrium distribution is again order of h'^ in the 
exponent. Writing the equilibrium distribution in the form 

P{x, t) = N exp { [-V{x) + h'^F{x)] / D} , (76) 

the function F[x) satisfies the equation [V = V{x), F = F{x) 

= SGD'-'F' + 6V'^ - GODV'^V" + 96£>W"2 + TSD^V^' V" 

-90D^V"V"' - 66D^V'V^^^^ + 2iD'^V^^^ . (77) 

It is clear that the two best algorithms, in the sense that they generate an 
equilibrium distribution which is the closest to the "true" one, are the Heun 
scheme and the "full algorithm" . It should be appreciated that even if the 
correction in the equilibrium distributions may look tiny (one could always 
think of making the integration time step smaller), the actual Fokker-Planck 
equation is different, and the difference in the equilibrium distribution is in an 
exponent. Hence, we may expect unpredictable results when things like the 
mean first passage time between two minima in a bistable potential are com- 
puted. What is worse, for fixed h the correction depends on D, making it more 
difficult to compute, for instance, activation energies from the simulations. 

To show that the calculation of the Mean First Passage Time (MFPT) is 
indeed strongly infiuenced by the algorithm chosen, I summarized the result of 
computing the mean first passage time between two minima in the potential 
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V{x) = + using different algorithms and integration time steps in 

figure 3. The noise intensity was always D = 0.1: note how the convergence to 
the "limiting" MFPT is much slower when I used algorithms which lead to an 
order h correction in the equilibrium distribution (a similar conclusion could 
have been guessed from figure 3). I used 10000 averages for each point in figure, 
so that the statistical scattering is order of 1.0%. Overall, the full algorithm is 
the one which stays closer to the limiting MFPT as the integration time step 
is changed, with the Heun algorithm which is typically also very close. The 
expected theoretical value for the mean first passage time is 66.3. 

I will finish this section with yet another example which should warn us 
about the care which must be exerted when higher order integration schemes 
are derived. One of the earliest, and most "transparent", derivation of in- 
tegration schemes, based on Runge-Kutta, is the one in [53]. There a first 
algorithm is derived as an example (termed 2o2s1g by the author, meaning 
a second order deterministic one, second order stochastic one and needing the 
generation of one gaussian deviate), and a second algorithm is then given in 
appendix (3o3s2g, supposedly a better and more accurate algorithm), as a 
"production" algorithm. 

Now, whereas the "example" algorithm of [53] is equivalent to the Heun 
scheme, and it shares all "good" properties of this algorithm, the 3o3s2g 
scheme leads to the equilibrium distribution 



i.e. it generates an equilibrium distribution which differs by the correct one 
by an order h in the exponent, and we should expect problems similar to the 
ones typically found in very low order algorithms. Also, it apparently does not 
do much better than the Z^-less scheme, which requires one less calculation of 
the deterministic force for each integration step. 

3 Higher order stochastic differential equations 

3.1 One pole and two poles filters 

An important role is played by noisy drivings which are not white, but with a 
spectral density which is Lorentzian. In particular, systems where the stochas- 
tic forcing is first passed through a one pole or a two poles filters. If f{t) is the 
usual white noise gaussian process, and y{tO is the noise driving the system of 
interest, y{t) is described by the equations 



P{x,t) = Nexp{[-V{x) - h {0.14:5186DV" {x})] / D} . 



(78) 



y = --y + f{t) 



(79) 
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y = -jy - ^ov + fit)- 



(80) 



The linearity of these equations has been exploited in [47] and [54,55] to derive 
an algorithm which parallels the algorithm of Eq. 12, with suitable Zi , and Z2. 
Also in this case we can improve the basic algorithm using the Heun algorithm 
(Eq. 39). The important feature of these algorithms is that they are fully 
implicit for the integration of the noise: this implies that they will exactly 
reproduce the evolution of the filtered noise, irrespectively of the relative ratio 
between the integration time step and the intrinsic time scales of the filters. 
This is confirmed using the algorithm for a one pole filter to integrate a quasi 
white noise: in figure 3 I plotted the result, and we should note that for all dh 
considered, dh is larger than r. 

3.2 Some alternative algorithms 

For the specific case of exponentially correlated noise it is possible to find in 
the literature some alternative integration scheme. Here I will briefiy review 
the one introduced in [7,9] and the one introduced by [10]. 

A very elegant and efficient scheme for exponentially correlated noise has 
been proposed in [7,9]. The idea is to start from 

x = f{x)+y (81) 
with y gaussian and correlated as per 

{y(t)y(s))=F(t-s). 

Clearly, we have through a Fourier transform 

{y{oj)y{oj')) = G{oj - oj') = H{oj)5{oj - oj'). 

Now, write (central limit theorem) 

N 

y{jh) = (l/Vh) «i cos{iAi^jh + <f>i) (82) 

i = 

with (j>i uniformly distributed over [0, 27r], a? = H{i Ah')Ah', and suitably cho- 
sen N and Ai>, such that the conditions 

Aiy < ly, 

and 



22 



are satisfied. It is tlien possible to integrate Eq. 81. The key point (see [7]) 
is that the sum appearing in Eq. 82 can be evaluated with a Fast Fourier 
Transform, which implies that to generate, say, y{ih) for times from up to 
Nh, only NlogN operations are required. Although for one pole or two poles 
filters the direct simulation would be faster, for computationally expensive 
H[uj) this method is much preferable [8,7]. For vector and parallel machines, 
furthermore, the whole algorithm can be made to run very efficiently. 

Another algorithm to integrate exponentially correlated noise is the one 
proposed in [10]. The idea is to replace the exponentially correlated noise by 
a superposition of independent random telegraph processes. Suppose we have 
N random telegraph processes {sj}, oscillating between the states +1 and — 1, 
and such that 

(s,) = (s,.(t)s,.(0)) = %e-2«l*l. (83) 
Introduce the variable Y, defined as 

1 ^ 

>^=2^«., (84) 

i = l 

which takes values between —N 12 to #/2 in steps of unity. If M = Y + #/2 
is the number of s,- = 1 in the sum, then the variable Y changes its state to 
Y ^ Y±l with a probability M/N for Y ^ Y-1 and l-M/N for Y ^ Y + 1. 

It is possible to work out the statistical properties of the variable Y: one 
finds [10] 

(Y) = (Y(t)Y(O)) = (#/4)e-2«l*l, (85) 

and, given that Y is the sum of independent processes, the quantity Y / N be- 
comes a gaussian in the limit of large N . We have clearly a possible representa- 
tion of an exponentially correlated noise. The apparent advantage with respect 
to other approaches is that, between different times at which the random tele- 
graph signal switches, the noise is "constant", and one is left to integrate a 
standard differential equation. The method picks the time sequence {tk} at 
which the variable Y switches using the prescription 

tk+i=tk + Til{Na) (86) 

where ?y is a random variable exponentially distributed with a cutoff equals to 
one: clearly, the larger N (to achieve a better gaussianicity), the smaller the 
time steps generated. 

Apart from the problems connected with the determination of the optimal 
N (which is a difficult question, in practice, and can only perhaps be inferred 
in an indirect way, for instance checking for consistency of the results as N 
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is changed), it is not obvious at all that the method should yield faster algo- 
rithms or, more important, that it reproduces distributions which are close to 
the correct one. In fact, taking a concrete examples, I integrated with this 
algorithm the SDE corresponding to the usual bistable potential, for r = 1.0 
and _D = 0.1. I followed the stochastic trajectory for a time equals to lO'*. The 



Table 1: Times taken by different algoritlims 



Time step 


Modified Heun 


Present Algorithm 


0.1 


11.17 


11.00 


0.05 


22.13 


14.53 



results are summarized in table 1, where I used # = 21 as suggested by the 
authors of [10]. At first sight, the present algorithms is better as I decrease 
the integration time step (although, as I have shown before, the Heun scheme 
works even for time steps as large as 0.1). The problem, however, is when I 
check the gaussianicity of the noise generated by the algorithm. Given the cor- 
relation time of the noise and the total integration time, I would expect that the 
statistical scattering on the second moment of the noise is around 1.0%. The 
result of calculating the various moments, for an integration time step equals 
to 0.05 (the most favorable case for the present algorithm), is summarized in 
table 2. It is clear that the present algorithm generates a noise which has a 



Table 2: Noise moments generated by the different algorithms 



Quantity 


Modified Heun 


Present Algorithm 


Theoretical 


(y) 


-2.9810-3 


-1.6310-3 


0.0 


{y') 


1.028 


1.265 


1.0 


{y') 


-7.6710-3 


-8.5110-2 


0.0 


{y')/({yY) 


3.063 


2.299 


3.0 


{y')/m{y')) 


5.224 


3.766 


5.0 



distribution which is very different from a gaussian one: the second moment is 
some 25% larger than expected, whereas the higher moments are much smaller 
than what they should be. In fact, a much more reasonable value for N , such 
that the noise has the correct distribution (on this time scale) is N ^ 2000, 
which would, however, make the whole algorithm some factor 100 slower. 
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3.3 Quasi conservative systems 

Another important class for which dedicated algorithms can be derived is given 
by the equations of motion 



ij = -jv + F{x)+£,{t) (87) 
where ^(t) is a random gaussian noise, with zero average and standard deviation 

mm) = ^iDS{t-s). (88) 

The above equation is commonly found in the liquid state literature, and 
several algorithms have been proposed, over the years, for its integration. The 
definitive word is perhaps summarized in [56], where two algorithms are pro- 
posed (see also references therein). The first algorithm integrates the above 
equation using the prescription 

x{t + h) = x{t) + cihv{t) + C2h'^F{x{t)) + t]i (89) 
v(t + h) = cov(t) + cihF(x(t)) + r]2, (90) 



where 



CO = e--"' (91) 

= ^ ('^^' 

and where rji and ri2 are two random gaussian variables with zero averages and 
moments 

(ryl) = D{1- e-'^'') (95) 

(,yi,y2) = -(l-e-^'')'. (96) 
7 

It is possible to check which is the equilibrium distribution reproduced by 
this algorithm. We know that the theoretical equilibrium distribution of Eq. 87 
should be given by 

P{x,v) = Nexp{-[v'^/2+V{x)] /D} (97) 
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where N is a normalization constant and V{x) = — J F{x)dx. Writing the 
distribution generated by the numerical scheme in the form 



P{x,v) = Nexp{- [vy2 + V{x) + hS{x,v)] /D} (98) 
we have that S{x,v) satisfies the difli"erential equation 

d^S{x,v) V dS{x,v) f F{x) _ v \ dS{x,v) ^ 1 ^ ^ 

(9t;2 D-) dx \ D-) Dj dv 27 

(99) 

This means that this algorithms fails to reproduce the correct equilibrium 
distribution at 0(h) in the exponent. It is also possible, for the case when 
V{x) = oj^x^ /2, to derive the equilibrium distribution at lower order in h, 
which reads 

P(x, v) = Nexp{- [v'-'/2 + w2a;2/2] /£)'} (100) 

with 

= (101) 

which means that when 7 goes to zero, the effective temperature simulated 
by the algorithm goes to zero (it must be said that in [56] it is acknowledged 
that the algorithm does not work well in this limit, although no formal proof 
is provided). 

To overcome the problems with the case of small 7, in [56] a second algo- 
rithm is proposed, which reads 

x{t + h) = x{t) + cihv{t) + C2h'^F{x{t)) + t]i (102) 
v(t + h) = cov(t) + (ci - C2)hF(x(t)) + C2hF(x(t + h)) + r]2, (103) 

which indeed reproduces the correct equilibrium distribution at o{h) (of course, 
there are corrections at 0{h^) in the equilibrium distribution, but I will not 
write them here). 

However, it is clear that both algorithms are fairly expensive in terms of 
CPU times: they require at least two random gaussian deviates and one or two 
evaluations of the force at each integration time step. Beside the Heun scheme 
and the full algorithm, which could be applied also to this case, it should be 
possible to find a more efficient algorithm for this specialized case. 

The basic idea is that if we took 7 = in Eq. 87, I could integrate this 
equation using a simplectic scheme (see, for instance [57]). The scheme which 
is straightforwardly used in the presence of the noise would integrate Eq. 87 
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as 



i = x(t) + ^v(t) (104) 

v(t + h) = v(t)+hF(£) (105) 

x(t + h) = x+^v(t + h). (106) 

It is then possible to reintroduce both the dissipation {—'jv) and the noise, to 
obtain the scheme 

i = x(t) + ^v(t) (107) 

v{t + h) = C2 [civ{t) + hF{i) + diri] (108) 

x{t + h) = x+y^v{t + h), (109) 

where ?y is a gaussian variable, with standard deviation one and average zero, 
and 

ci = l-f (110) 

C2 = (111) 



di = ^/2Djh. (112) 

This last algorithm does reproduce the correct equilibrium distribution 
(apart from terms o[h'^), as the second algorithm proposed in [56]), but it only 
requires one gaussian deviate and one evaluation of the force. Furthermore, 
in the limit 7 — 0, by construction it conserves the energy of the system, at 
o{h^). 

4 Conclusions 

I have shown algorithms to integrate general stochastic differential equations, 
discussing connected problems. I would like to add that any simulation in 
stochastic dynamics requires a good and fast random number generator. Very 
many generators are available in the literature: personally, after trying a few, 
I settled for a couple of algorithms which satisfied my needs, and which have 
the good features that are fairly fast and portable to all platforms (they are 
written in high level languages). For the generation of fiat distributions, I use 
the subtract and borrow algorithm, first proposed by Marsaglia (RCARRY), 
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in the implementation by Luscher [58,59]. For the generation of gaussian or 
exponential deviates, I use the Ziggurath algorithm, by Marsaglia and Tsang 
[60], modified to use the output from the subtract and borrow routines. Let me 
finish noticing that no simulation can be better than the random noise used: 
so, extreme care should be given to the choice of the random noise generator. 
Most library random noise generators are reasonable for small simulations, but 
when a large number of random terms is needed, or the gaussianicity and true 
randomicity of the noise is required, my strong advice is to test very carefully 
the algorithm one is going to use. To this end, the benchmarks proposed in [61] 
are still the state of art. 
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